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Computer simulations of first-order phase transitions using "standard" toroidal boundary condi- 
tions are generally hampered by exponential slowing down. This is partly due to interface formation, 
and partly due to shape transitions. The latter occur when droplets become large such that they 
self-interact through the periodic boundaries. On a spherical simulation topology, however, shape 
transitions are absent. By using an appropriate bias function, we expect that exponential slow- 
ing down can be largely eliminated. In this work, these ideas are applied to the two-dimensional 
Widom-Rowlinson mixture confined to the surface of a sphere. Indeed, on the sphere, we find that 
the number of Monte Carlo steps needed to sample a first-order phase transition does not increase 
exponentially with system size, but rather as a power law r oc V a , with a w 2.5, and V the system 
area. This is remarkably close to a random walk for which orw = 2. The benefit of this improved 
scaling behavior for biased sampling methods, such as the Wang-Landau algorithm, is investigated 
in detail. 
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I. INTRODUCTION 

Phase transitions in colloidal suspensions are of pro- 
found practical interest. Think, for instance, of phase 
separation in colloid-polymer mixtures [T], or the freez- 
ing of colloidal hard spheres at high densities [5] . For this 
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FIG. 1: Phase separation snapshots of the 2D WR mix- 
ture as obtained in grand canonical MC simulations using 
"standard" toroidal boundary conditions, i.e. a square with 
periodic boundaries. The snapshots were obtained at fugacity 
z = 2.5, which is well above the critical fugacity z CI , and so the 
transition is strongly first-order here. The light regions cor- 
respond to the vapor phase, dark regions to the liquid phase. 
Starting in the vapor phase (a), a droplet of liquid nucleates 

(b) . The droplet grows until the strip configuration is reached 

(c) . Further increasing the amount of liquid phase leads to a 
droplet of vapor (d), and finally this droplet vanishes, leading 
to the pure liquid phase (e). 



reason, there is also an enormous interest in the mod- 
eling of phase transitions by means of computer simula- 
tion. The investigation of phase transitions via computer 
simulation is not trivial, as there are numerous hurdles 
to overcome. One obvious problem is the issue of finite 
system size. Since computational resources are limited, 
one always deals with finite numbers of particles, whereas 
phase transitions are defined in the thermodynamic limit. 
Hence, there is an obvious "gap" to bridge, achieved in 
practice using finite-size scaling 

Another problem, which we focus on in the present pa- 
per, concerns exponential slowing down, and affects sim- 
ulations of first-order phase transitions [I]. To illustrate 
the problem we consider phase separation in the Widom- 
Rowlinson (WR) mixture [5] in two dimensions (2D). In 
this model there are two particle species, A and B, each 
modeled as disks of diameter a (in what follows a will be 
the unit of length). The only interaction is a hard-core 
repulsion between A and B disks. As is well known, the 
WR mixture phase separates provided the fugacities za 
and Zb, of A and B particles, are high enough; due to 
symmetry it holds that za = zb = z at the transition, 
which we shall use throughout this work. When phase 
separation occurs, one obtains a "vapor" phase (dense in 
B species, and lean in A species) and a "liquid" phase 
(lean in B species, and dense in A species). The particle 
density of A species pa — Na/V may be used as order 
parameter to distinguish between the phases, since it as- 
sumes a low value in the vapor phase, and a high value 
in the liquid. Here, Na is the number of A particles in 
the system, and V the total system area (the choice for 
A or B is arbitrary, of course). Despite its simplicity, the 
2D WR mixture is relevant for phase separation in cell 
membranes, as the latter also constitute effectively 2D 
systems (which even give rise to 2D Ising critical expo- 
nents [5]). 
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The standard approach to simulate the 2D WR mix- 
ture is to perform a grand canonical Monte Carlo (MC) 
simulation on a V = L x L square with periodic bound- 
aries. Provided the fugacity z significantly exceeds the 
critical value z cr , such that the transition is first-order, 
phase separation proceeds as shown in Fig. [T] [7j [8] . Start- 
ing in the vapor phase (a), a nucleation event occurs, 
leading to the condensation of a droplet of the liquid 
phase (b). The droplet grows until it interacts with it- 
self through the periodic boundaries, leading to the strip 
configuration (c). In the strip configuration, vapor and 
liquid coexist with each other, separated by two inter- 
faces that run perpendicular to one of the edges of the 
simulation square as this minimizes the intcrfacial area. 
The approach to the liquid proceeds via the formation of 
a vapor droplet (d), which eventually vanishes, leading 
to a pure liquid phase (e). 

The path connecting vapor and liquid thus passes the 
strip configuration of Fig. [TJc) . However, in a stan- 
dard MC simulation, where configurations appear pro- 
portional to their Boltzmann weight, the strip configura- 
tion is extremely rare, due to the large amount of inter- 
face that it contains. In 2D, the total interface length in 
the strip configuration equals 2L, corresponding to a free 
energy barrier AF = 2aL, where a is the line tension. 
Hence, starting in one of the pure phase (a) or (e), it typ- 
ically takes r ~ exp (2a L) MC steps to reach the strip 
configuration. Since the "tunneling" time r increases ex- 
ponentially with the system size L, and hence explains 
the phrase "exponential slowing down" , it is clear that 
the standard MC method must be modified at a first- 
order phase transition in order to remain efficient. 

Such modifications have been made, and the state-of- 
the-art is to not sample from the Boltzmann distribution, 
but from a modified distribution, such that the "unfavor- 
able" interface configurations of Fig.flTb,c,d) are sampled 
just as often as the "pure" phases (a) and (e). Crucial 
to these methods is the use of an order parameter, con- 
structed such that 

(1) it varies strictly monotonically along the path con- 
necting the phases, and 

(2) is computationally fast to calculate. 

Regarding the 2D WR mixture, a convenient order pa- 
rameter is the particle density pa, which certainly fulfills 
criterion (2). In the vapor phase pa = PA.vap, in the 
liquid phase pa = /?A,iiq, while in the strip configuration 
Pa ~ (PA,va P + PA,iiq)/2- The idea is to perform a MC 
simulation using a biased energy Eb = Eq + w(pa), with 
Eq the original energy of the system, and w(pa) some 
a priori unknown function of the order parameter pa- 
Clearly, by tuning w(pa) appropriately, the probability of 
the interface configurations can be artificially enhanced. 
The aim is to construct w(pa) such that the simulation 
performs a random walk in pa- Various methods can be 
used to construct w(pa) in practice, such as multicanon- 
ical sampling jl] , (successive) umbrella sampling [TU] , 
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FIG. 2: The analogue of Fig. [T]but this time the simulation 
was performed on the surface of a sphere. In contrast to 
toroidal boundary conditions, the transition from the vapor 
(a) to the liquid (e) involves only the nucleation of a droplet. 
Once such a droplet has formed, the path from (b) to (d) does 
not involve any shape transitions; the images (b), (c), and (d) 
only differ in the relative volume of the phases, not in their 
geometric configuration. 



and Wang-Landau sampling [11 . The methods differ in 
details, but all provide a means to obtain w(pa)- 

Hence, using a suitable w(pa), the free energy bar- 
rier of interface formation is eliminated, the coexistence 
configurations of Fig. [ljb,c,d) become accessible, and a 
random walk in pa will result (or so one hopes). In prac- 
tice, this is not the case because, strictly speaking, pa 
does not fulfill criterion (1) and hence is not a suitable 
order parameter. To see this, consider the case where 
half of the simulation square is occupied with vapor, and 
the other half with liquid. One way to arrange the phases 
is the strip configuration (c), with the distance between 
the interfaces being L/2. However, one could equally well 
arrange the phases in one of droplet configurations (b) or 
(d), with the droplet radius being L/^/2tt. Both "solu- 
tions" yield the same order parameter pa, but clearly dif- 
fer in topology. This is a problem because the transition 
from the droplet to the strip configuration is a first-order 
transition by itself, with a complicated order parame- 
ter not simply related to pa 0- As with any first-order 
transition, these so-called shape transitions also lead to 
exponential slowing down. This means that simulations 
do not yield random walk behavior in the order parame- 
ter, even when w(pa) is accurately known. Instead, most 
time is spend in the pure phases (a) and (e), or in the 
strip configuration (c), but transitions between the pure 
phases and the strip become increasingly rare with in- 
creasing system size 0112]. This leads to poor sampling 
statistics in practice. 

In principle, these problems can probably be overcome 
using a more sophisticated order parameter, one that also 
distinguishes between shape. However, such order pa- 
rameters are not trivial to construct, and are likely to 
be computationally expensive, which would violate cri- 
terion (2). An alternative approach is to use a different 
simulation box topology, one in which shape transitions 
are absent . Note that the droplet-strip transition 
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is a consequence of using a square simulation box with 
periodic boundaries (topologically equivalent to a torus) . 
It is clear that on the surface of a sphere, the droplet- 
strip transition would not occur. Instead, on a sphere, 
we expect a first-order transition to proceed as shown in 
Fig. [2] Starting in the vapor phase (a), we still expect a 
nucleation event, leading to a droplet of the liquid phase 
(b). By increasing pa further, the droplet grows contin- 
uously (c,d) but there are no shape transitions. There is 
still one nucleation event, of course, as the vapor droplet 
(d) vanishes, leading to a pure liquid phase (e). Hence, 
on the surface of a sphere, shape transitions are absent, 
and by using an accurate bias function w(pa) one should 
achieve behavior more closely resembling a random walk 
in pa- In fact, any remaining slowing-down is due to nu- 
cleation, i.e. of actual physics taking place, and should 
prove rewarding to study further. 

The primary aim of this paper is to investigate if expo- 
nential slowing-down at first-order transitions is indeed 
eliminated on a spherical topology. The idea of doing so 
was announced some time ago [7], but as far as we know, 
such simulations have not been performed to date. In 
fact, we are only aware of Ref. which considers a 2D 
Ising model on the surface of a cube. The latter only 
crudely approximates a sphere, but already the barriers 
arising from shape transitions were seen to soften. Of 
course, being a lattice model, the extension to a spheri- 
cal topology is not feasible in the Ising case. In contrast, 
for the 2D WR mixture, which is entirely off-lattice, the 
extension to a spherical topology poses no fundamental 
objections. 

The outline of this paper is as follows. We first de- 
scribe the grand canonical MC method in the presence 
of a bias function, and we provide some implementation 
details as to how this method can be efficiently imple- 
mented on the surface of a sphere. Next, we perform a 
number of cross-checks, to demonstrate that toroidal and 
spherical simulation topologies are consistent with each 
other in the thermodynamic limit. We then turn to our 
main result, and show that at first-order phase transi- 
tions, exponential slowing down on the sphere is largely 
eliminated. Finally, we discuss how this improved per- 
formance benefits a number of sampling algorithms, in 
particular the Wang-Landau algorithm. We end with 
some concluding remarks in the last Section. 



II. METHODS 
A. Grand canonical Monte Carlo 

We simulate the 2D WR mixture in the grand canoni- 
cal (GC) ensemble using a bias function w(Na), i.e. de- 
fined on the number of A particles. The WR mixture 
was defined in the Introduction; here we only explain the 
GC simulation method. In the GC ensemble, the fugac- 
ity z and the system area V are fixed, while the number 
of particles fluctuates. Phase separation is studied using 



the order parameter distribution Pv(Na\z), defined as 
the probability to observe a system containing Na par- 
ticles of type A. We emphasize that Pv(Na\z) depends 
on the imposed fugacity z, the system area V, and on 
the simulation box topology (here: toroidal or spheri- 
cal) . The basic MC steps used to sample the distribution 
are insertions and removals of single particles. At each 
step, the simulation attempts with equal probability one 
of the four following moves: 

(1) insertion of one A particle at a random location, 

(2) removal of one randomly selected A particle, 

(3) insertion of one B particle at a random location, 

(4) removal of one randomly selected B particle. 

If the insertion attempts lead to forbidden overlaps, they 
are rejected. Otherwise, the moves are accepted with 
probabilities 
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In the above, Na (Nb) refers to the number of A (B) 
particles in the system at the start of the move. Note the 
presence of the bias function w(Na) in moves involving A 
particles. As was explained in the Introduction, the bias 
function is needed to overcome the free energy barrier of 
interface formation. 



B. Implementation 

The most CPU consuming steps are particle insertions, 
since here one needs to check for overlap with particles of 
the opposite species. We now discuss how these checks 
can be performed efficiently on the surface of a sphere. To 
simulate a total area V, the sphere radius must be R = 
v/V '/47r. The position of each particle on the sphere is 
stored using a 3D vector r = (x, y, z) with \ f\ = R. This 
means carrying around a third coordinate but allows to 
eliminate time-consuming trigonometric functions. First 
note that the on-sphere distance d between two particles 
i and j is the length of the shortest path over the sphere. 
This path lies on a great-circle, i.e. a circle with radius 
R. Hence, d = R9, where cos 9 = {fi ■ rj) / R 2 . Unlike 
particles i and j overlap when d < a, with a the particle 
diameter or, equivalently, whenever 



> R 2 cos (a/R) 



(5) 
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FIG. 3: Order parameter distribution Pv(Na\z) for the 2D 
WR mixture obtained at fugacity z — 2.5 using toroidal 
(dashed line) and spherical (solid line) simulation topologies. 
Note the logarithmic scale. The system area equals V — 1600 
in both cases. Due to the high value of z, the vapor peak on 
the left is squeezed to the edge. Also indicated are the free 
energy barriers AF t and AF 3 , which can be used to obtain 
the line tension a. 



where the right-hand side is a constant, which needs to be 
evaluated only once at the start of the simulation. Hence, 
to check for overlap, only the computationally cheap term 
r*j • fj is needed but no trigonometric functions. 

The second optimization concerns the implementation 
of link-cell neighbor lists [T3] on the sphere. To this end, 
the sphere (of radius R) is "embedded" in a 3D cube 
of edge 2R. The cube itself is partitioned iufixnxn 
equally sized sub-cubes, with n = 2R/a rounded down 
to the nearest integer. Since it holds that 

d>\n-r,l (6) 

with d the on- sphere distance between particles i and j, 
one only needs to check for overlap with particles that 
are in the same sub-cube, or in any of the neighbor- 
ing sub-cubes (including diagonal neighbors). Note that 
the number of neighboring sub-cubes that needs to be 
checked is typically less than the maximum of 3 3 — 1 = 26 
possible neighbors, since only sub-cubes actually inter- 
secting with the surface of the sphere have to be taken 
into account. In practice, only about 13— 14 neighboring 
sub-cubes were counted in our simulations. 



III. RESULTS 

A. Order parameter distribution and line tension 

We begin our analysis by explicitly showing the or- 
der parameter distribution Pv{Na\z) as obtained using 
a toroidal and spherical simulation topology. Fig.[3]shows 
a typical result for a fugacity z high above the criti- 
cal fugacity z cr ; the system area V is the same in both 
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FIG. 4: Variation of the line tension a versus z as obtained 
using toroidal (dashed line with squares) and spherical (solid 
line with circles) simulation topologies. The data were ob- 
tained at fixed system area V = 1600 in both cases, and so 
finite-size effects are not accounted for. 



cases. The distributions reveal two peaks: the left peak 
corresponds to the pure vapor phase, the right peak to 
the pure liquid, and from the peak positions pA,vap and 
PA,iiq can be read-off. Whenever the simulation visits the 
peaks, a single homogeneous phase is observed, i.e. re- 
sembling the snapshots of Fig.[lja,e) and Fig.[2ja,e). On 
the scale of the graph, the peak positions between the 
toroidal and spherical topology practically coincide. This 
is to be expected because the pure phases do not contain 
any interfaces. Only in the region between the peaks 
are differences between the two topologies expected to 
appear. 

For the toroidal topology, a pronounced flat region be- 
tween the peaks unfolds. This is where the system as- 
sumes the strip configuration of Fig. [ljc). In this config- 
uration, an increase of the volume of either phase at the 
expense of the other phase merely moves the interface 
but does not change its form nor affect the bulk phases. 
Hence, the free energy remains the same under such a 
change which is the origin of the characteristic flat region 
in the probability distributions of systems on a toroidal 
topology. Following Binder [T3], the average height AF t 
of the peaks above the flat region, measured in the log- 
arithm of Pv(Na\z), corresponds to the free energy cost 
of interface formation. Since the total amount of inter- 
face in the strip configuration equals 2L, one immediately 
obtains the line tension 

a = AF t /(2L) (toroidal topology), (7) 

with L the edge of the simulation square. In contrast, 
using a spherical topology, Pv(Na\z) does not reveal any 
flat region. The reason is that on the sphere, any change 
in the relative volume of the phases inevitably creates or 
destroys interface. The maximum amount of interface is 
generated when half the sphere is occupied with vapor, 
and the other half with liquid, i.e. conform Fig. [2lc). 
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FIG. 5: Determination of z cr via finite-size scaling using both 
toroidal (squares) and spherical (circles) simulation topolo- 
gies. Plotted are the fugacities at which x an d Y^ 1 attain 
their extrema versus l/l, with I = yfV the lateral extension of 
the system; the lines are linear fits. For systems of finite size, 
the results differ significantly between the different topologies, 
but agree on the value of z CI in the thermodynamic limit. 



The interface length then equals 2ttR, and the analogue 
of Eq.([7]) becomes 



a = AF s /(2irR) (spherical topology), 



(8) 



with R the sphere radius, and AF S the barrier height. 
The estimates of the line tension as a function of the 
fugacity are shown in Fig. [4] for both topologies using 
V = 1600. In agreement with theoretical expectations, a 
decreases with decreasing z, and at z cr it vanishes. How- 
ever, <7 obtained on the sphere is systematically below 
that of the torus, the discrepancy being around 5%. In 
principle, we only expect agreement in the thermody- 
namic limit, and so a detailed investigation of finite-size 
effects [15] is required to resolve this issue. In addition, 
a obtained on the torus corresponds to planar interfaces, 
whereas the interfaces on the sphere are curved. Hence, 
there could be curvature corrections, possibly involving 
Tolman's length [To] . 



B. Locating the critical point 

In the thermodynamic limit, phase transition proper- 
ties should not depend on the simulation topology. So 
as a test for the validity of using a spherical topology, 
rather than the more common toroidal one, we compare 
predictions for the critical point. Simulating at fugacities 
around the critical point, and using histogram reweight- 
ing [T7], we verified that in both cases the same value 
of the critical fugacity z cr is obtained, as well as critical 
exponents consistent with 2D Ising universality. Fig. [5] 
shows the result of a finite-size scaling study along the 
lines of Ref. UHl For both topologies, we have plotted the 
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FIG. 6: Main result of this paper: shown is the number r 
of MC steps needed to traverse from pA — to pa = 3.5 
and back as a function of the system area V in biased simu- 
lations using toroidal (dashed line) and spherical (solid line) 
boundary conditions. Note the double-logarithmic scale! The 
important result to take from this figure is that r increases 
strongly with V (presumably exponentially) on the torus, but 
only weakly (power law) on the sphere. The data were ob- 
tained for the 2D WR mixture at fugacity z — 2.5; the ex- 
ponent of the power law for the spherical topology equals 
a w 2.5, which is remarkably close to «rw = 2 of a true 
random walk. 



fugacities at which the susceptibility \ attains its maxi- 
mum versus l/l. Here, I = \fv denotes the "length" of 
the system. Also shown are the fugacities at which the 
generalized susceptibility attains its global extrema, 
with Y^ defined in Ref. \M The lines arc linear fits, 
which capture the data well, consistent with the exact 
2D Ising value v = 1 of the correlation length critical 
exponent. Clearly, in finite systems, the results between 
the two topologies differ, but both extrapolate to a com- 
mon value z cr = 1.717(2) in the thermodynamic limit 
(the error reflects the scatter between individual scaling 
results). Furthermore, in both topologies, we find that 
the maximum value of the susceptibility increases with 
the length of the system as Xmax <x P^", with 7 the crit- 
ical exponent of the susceptibility. We obtain j t w 1.754 
and 7 S w 1.743, for the toroidal and spherical topology, 
respectively, which both compare well to the exact 2D 
Ising value 7 = 7/4. Our estimate of z CT for the 2D WR 
mixture is close to the one reported in Ref. 1191 but upon 
careful inspection does underestimate it (by about 0.5%). 
Interestingly, one of us (RV) has experienced a similar 
disagreement for the 3D WR mixture as well 20J . The 
origin of the discrepancy is not clear. 



C. The fate of exponential slowing down 

We now come to the main result of this paper, where we 
consider exponential slowing down at a first-order phase 
transition. Our hope is that, by using a spherical topol- 
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TABLE I: Relative accuracy of selected physical observables 
obtained using WL sampling on spherical and toroidal topolo- 
gies for several system areas V. The ratios all exceed unity, 
implying that the spherical topology is to be preferred. The 
data were obtained for the 2D WR mixture at fugacity 
z = 2.5. 



ogy, exponential slowing down can be largely eliminated. 
To this end, we set the fugacity to z — 2.5 which is well 
above z CT , and so the transition is strongly first-order. 
We remind the reader that our simulations use a bias 
function w(Na) to overcome the free energy barrier of 
interface formation, and also that w(Na) is a priori un- 
known. Hence, for a number of system sizes V, accurate 
bias functions w(Na) were first obtained using successive 
umbrella sampling [10J, for both toroidal and spherical 
topologies. Next, biased simulations were performed us- 
ing the (now known) bias functions, and the number of 
MC steps t needed to traverse from pa = to pA = 3.5 
and back was measured; the reader can verify in Fig. [3] 
that this range is sufficient to sample both the vapor and 
liquid peaks. 

The resulting r data are collected in Fig. [6] The most 
striking feature of the plot is that for a toroidal topology 
r indeed increases faster than a power law with V . Fol- 
lowing the discussion in the Introduction, we attribute 
this slow down to free energy barriers associated with 
the droplet-strip shape transition, which are not over- 
come by the bias function w(Na)- The data for the 
spherical topology, in contrast, do not reveal any expo- 
nential slow down for the system sizes considered here 
but rather a power law increase; approximately r ~ V a , 
with a ~ 2.5. This is still a slow down compared to a 
perfect random walk, for which a^w = 2, but not an ex- 
ponential one. Possibly, the remaining slow down is due 
to nucleation events. Comparing the values of r between 
the two topologies, it is striking that already for system 
area V = 1600, r on a torus is ten times that of r on a 
sphere. 
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FIG. 7: Number of MC moves required to complete one WL 
iteration, using a toroidal (vertical bars) and spherical (cir- 
cles) simulation topology, on double logarithmic scales. Re- 
sults are shown for system areas V = 900, 1600, 2500 (bottom 
to top). Note that during WL sampling the larger values of 
5 are sampled first. The data were obtained for the 2D WR 
mixture at fugacity z = 2.5. 




FIG. 8: Variation of the statistical error E in Aw(Na) versus 
Pa, as obtained in WL simulations on toroidal (dashed line) 
and spherical (solid line) topologies. For clarity, intervals of 
20 error estimates are combined to a single average. The 
error bars on the toroidal data represent jackknife errors for 
the average of this re-binning; error bars for the spherical 
topology are omitted. The data were obtained for the 2D 
WR mixture at fugacity z — 2.5 and system area V = 2500. 



D. Wang-Landau sampling 

Having shown that the tunneling time r using a spher- 
ical topology can be much smaller compared to that on 
a torus, a natural next question is whether this improve- 
ment also increases the performance of the algorithms 
used to construct w(Na)- One such algorithm is Wang- 
Landau (WL) sampling 11J. Here, the bias function is 
initially set to w(Na) — 0, and one proceeds to simulate 
as explained in Section |II A| Each time a state with Na 
particles of type A is visited, the corresponding value of 
the bias function is decreased by a modification term: 



W(N A ) -> W(N A ) - S. This reflects the idea that the 
current number of A particles is just being found a bit 
more probable than previously expected and hence needs 
a smaller bias. One continues to simulate until all par- 
ticle numbers Na over the range of interest have been 
visited sufficiently often [5T] , which completes one WL it- 
eration. At this point, the modification term is reduced, 
say, 5 — > (5/2, and the next iteration is started. By using 
a large modification term at first, say 5=1, one ensures 
that all states will be visited in relatively short times. At 
later stages, when 6 is small, changes to the bias function 
become negligible, and the algorithm is said to have con- 
verged. Clearly, for the performance of this algorithm, 
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it helps if the simulation traverses the density range of 
interest as quickly as possible. This is obviously related 
to the tunneling time r, which is significantly reduced 
on a spherical topology, and so we expect an increased 
performance for WL sampling too. 

To test the performance of WL sampling on toroidal 
and spherical topologies, the bias function w(Na) was 
measured independently 40 times using system sizes V = 
900, 1600, 2500 (by independent we mean that each WL 
run was started with its own unique stream of random 
numbers). In Fig. [7] the average number of MC steps 
needed to complete one WL iteration is plotted as a func- 
tion of the modification term S. In the early stages, where 
5 is large, the number of MC steps is effectively identi- 
cal. In this regime S is so large that the simulations on 
the torus are still easily pushed through the regime of 
the droplet-strip transition, and hence there is no no- 
ticeable difference. However, at later stages, where S is 
small, the number of required MC steps is significantly 
less on the sphere, as expected, since here the sphere 
simulations benefit from the improved diffusion behavior 
demonstrated in Fig. [6] 

Hence, late-stage WL iterations indeed complete faster 
on a spherical topology, compared to a toroidal one. 
Next, it remains to be shown that actual physical ob- 
servables are also more accurately obtained. To this end, 
we consider the average density of A particles (pa), and 
the first four central moments 

m n = ( \p A - (pA)\ n ) ■ (9) 

Note that (pa) and m n are trivially obtained from the or- 
der parameter distribution Py(Na\z), which in turn is re- 
lated to the bias function w(N a) = — log Pv(Na\z). For 
each topology, we thus have a set of 40 estimates per ob- 
servable. Using the jackknife method [221, one can derive 
the statistical error S in each observable, and the "best" 
topology is the one with the smallest E. In Table [TJ the 
ratio of errors S t /E s is shown for each observable, with 
S 4 the error as obtained on the torus (sphere). The 
point to take from the table is that the ratios all exceed 
unity, meaning that the data from the spherical topol- 
ogy are more reliable. In combination with the findings 
of Fig. [7] we conclude that WL simulations on spherical 
topologies are overall more efficient. 

Finally, we demonstrate that the enhanced perfor- 
mance of WL sampling on spherical topologies can in- 
deed be attributed to the absence of the droplet-strip 
shape transition. To this end, we consider the difference 

Aw(N A ) = w(N A + 1) - w(N A ) (10) 

between adjacent weights. Again using the jackknife 
method, we estimated the statistical error S in Aw(Na) 
from the set of 40 simulations for each topology. Plot- 
ted in Fig. [8] is £ versus pa- The striking feature is that 
on the torus, £ displays two extra large peaks, which 
are completely absent on the sphere. These extra peaks 
in the torus data reflect the sampling difficulties aris- 
ing from the droplet-strip transition. Away from the 



V (pa) rxi\ m2 rm rrn 

900 1.61 2.61 2.62 6.8 1.83 

1600 3.47 6.34 6.34 14.6 2.67 

2500 4.43 11 11.1 30.1 11.2 

TABLE II: The analogue of Table [I] but this time using suc- 
cessive umbrella sampling [10] , 

droplet-strip transitions, both topologies yield essentially 
the same statistical error, as expected. Note also the ex- 
cellent agreement between both topologies regarding the 
peak on the far right of the graph: this peak reflects a 
sampling problem arising from nucleation, which indeed 
should occur in both topologies. In principle, a nucle- 
ation peak should also be visible on the far left of the 
graph, but due to the high fugacity used, this peak is 
probably squeezed onto the vertical axes. 



E. Successive umbrella sampling 

Another algorithm that can be used to construct the 
bias function w(Na) is successive umbrella sampling 
(SUS) [10, ■ When SUS is implemented on a spherical 
topology we also observe an increase in performance. In 
Table [H] we display some typical benchmarks regarding 
the accuracy of a number of physical observables. Com- 
pared to WL simulations, the performance increase in 
SUS simulations appears to be even more pronounced. A 
possible explanation may be that, unlike in WL sampling, 
each state Na was simulated using the same number of 
MC steps for both geometries. 

IV. CONCLUSIONS 

In conclusion, we have shown that the droplet-strip 
transition, which is inherent to toroidal systems under- 
going a first-order phase transition, acts as barrier mak- 
ing simulations of large systems increasingly harder. For 
the 2D WR mixture investigated here, the droplet-strip 
transition can be easily, and at the cost of a constant 
fraction of CPU time, be eliminated by simulating the 
system on the surface of a sphere. We have shown that 
for WL sampling [TT] and SUS [TO], simulations on the 
surface of a sphere yield better results. The droplet-strip 
transition is a general feature, also in 3D, appearing at all 
first-order phase transitions studied on toroidal topolo- 
gies. Hence, the use of a spherical topology is expected 
to be beneficial in a great number of other systems also. 

The 2D implementation sketched here should quite 
straightforwardly extend to 3D [53], and to other 
(short-ranged) interactions also; obvious candidates are 
the hard-core square-well fluid, the (cut-and-shiftcd) 
Lennard- Jones fluid, and colloid-polymer mixtures. 
However, models in which the pair potential is a more 
complicated function of the distance may require the use 
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of computationally expensive trigonometric functions, 
which were successfully circumvented in the present im- 
plementation. Even so, this additional computational 
effort should be outbalanced by the elimination of ex- 
ponential slowing down, provided the systems are large 
enough. 

Some models may not be so easily transferable to the 
sphere, in particular when particle orientation comes into 
play. An example is the 2D Zwanzig model [H] of hor- 
izontally or vertically aligned hard rods. While on the 
torus one can uniquely speak of horizontal and vertical 
directions, this is prevented by the intrinsic curvature 
on the sphere. Note also that the advantage of using 
a spherical topology is to eliminate exponential slowing 
down at first-order transitions. Around the critical point, 
where the transition is continuous, we did not see much 
advantage using the spherical topology. 

A different application where the use of a spherical 



topology may be beneficial is in the simulation of droplets 
|25j |26j [27] . One problem of using a toroidal topology 
is that the maximum size of the droplet that one can 
simulate is limited to the point where the droplet-strip 
transition takes place [37]. On a spherical topology, this 
problem is circumvented. Finally, we would like to point 
out that phase separation on the surface of a sphere is 
also realized experimentally in giant vesicles |28j ; confo- 
cal microscopy images of the latter qualitatively resemble 
the simulation snapshots of Fig. [2j 
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